accu = cell(numel(traindata), numel(regvals));
etimes = cell(size(accu));
fvals = cell(size(accu));
postdim = cell(size(accu));
macc = cell(size(accu));
mfval = cell(size(accu));
aucs = cell(size(accu));
mauc = cell(size(accu));

if any(~ismember({'rid_s','rid_t'}, who))
    rid_s = 1;
    rid_t = numel(regvals);
end

for di = 1:numel(traindata)
    fprintf('di = %d\n', di);

    y = traindata{di}.y;
    X = traindata{di}.X;
    tnidx = traindata{di}.nidx;
    if ~exist('validid', 'var') || ~validid
        ty = testdata{di}.y;
        tX = testdata{di}.X;
        ttidx = testdata{di}.nidx;
    else
        ty = validata{di}.y;
        tX = validata{di}.X;
        ttidx = validata{di}.nidx;
    end
    assert(tnidx(end)==numel(y)+1 && ttidx(end)==numel(ty)+1);
    
%     stdx = std(X(1:end-1,:),0,2);
%     scalex = 1./sqrt(invsigmawsqr)./stdx;
%     meanx = mean(X(1:end-1,:),2);
%     X(1:end-1,:) = (X(1:end-1,:)-meanx(:,ones(Ntrain,1))).*scalex(:,ones(Ntrain,1));
%     tX(1:end-1,:) = (tX(1:end-1,:)-meanx(:,ones(Ntest,1))).*scalex(:,ones(Ntest,1));
    
    for ci = rid_s:rid_t
        C = regvals(ci);
        fprintf('C = %.4f\n', C);
        C = C(ones(size(y))) + ((cratio-1)*C)*(y==1);
        
        Z = Z_init;
        W = W_init;
        eta = eta_init;
        invlambda = invlambda_init;
        
        if algtype ~= 2
            mzZ = zeros(D,0);
            mzW = zeros(Nall,0);
            mzeta = zeros(M,0);
            minvlambda = zeros(size(invlambda));
        end
        mscore = zeros(Ntest, 1);
        
        maxacc = 0;
        minfval = inf;
        loopi = numel(postdim{di,ci}) + 1;
        cnvg = 0;
        earlystop = 0;
        cmacc = zeros(2,1);
        cmauc = zeros(2,1);
        
        while true
            [Z,W,eta,invlambda,zcidx,newK,etime,fval] = mcmc(Z, W, eta, invlambda, ...
                [X,tX], y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, tnidx, ...
                C, ell, poisstrunc, algtype, loopi);
            K = size(Z, 2);
            
            pscore = zeros(Ntest, 1);
            for m = 1:M
                nind = ttidx(m):ttidx(m+1)-1;
                pscore(nind) = eta(m,:)*Z'*tX(:,nind);
            end
            cacc = mean(2*(pscore>0)-1 == ty);
            [~,~,~,cauc] = perfcurve(ty, pscore, 1);
            
            if algtype ~= 2 && loopi > burnin
                navg = loopi - burnin;
                if navg == 1
                    mnZ = Z;
                    mnW = W;
                    mneta = eta;
                else
                    [mnZ, mzZ] = updatem(mnZ, Z, navg, zcidx, newK, mzZ);
                    [mnW, mzW] = updatem(mnW, W, navg, zcidx, newK, mzW);
                    [mneta, mzeta] = updatem(mneta, eta, navg, zcidx, newK, mzeta);
                end
                mZ = [mnZ,mzZ]; mW = [mnW,mzW]; meta = [mneta,mzeta];
                minvlambda = minvlambda + (invlambda-minvlambda)./navg;
                
                mscore = mscore + (pscore-mscore)./navg;
                cmacc(2) = mean(2*(mscore>0)-1 == ty);
                [~,~,~,cmauc(2)] = perfcurve(ty, mscore, 1);
                
                pscore = zeros(Ntest, 1);
                for m = 1:M
                    nind = ttidx(m):ttidx(m+1)-1;
                    pscore(nind) = meta(m,:)*mZ'*tX(:,nind);
                end
                cmacc(1) = mean(2*(pscore>0)-1 == ty);
                [~,~,~,cmauc(1)] = perfcurve(ty, pscore, 1);
                
                cmfval = fobj(mnZ, mnW, mneta, [X,tX], y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, tnidx, C, ell);
            else
                mZ = Z; mW = W; meta = eta; minvlambda = invlambda;
                if algtype == 2
                    cmacc(1) = cacc; cmfval = fval; cmauc(1) = cauc;
                    mscore = mscore + (pscore-mscore)./loopi;
                    cmacc(2) = mean(2*(mscore>0)-1 == ty);
                    [~,~,~,cmauc(2)] = perfcurve(ty, mscore, 1);
                else
                    cmacc(:) = nan; cmfval = nan; cmauc(:) = nan;
                end
            end
            
            accu{di, ci} = [accu{di, ci}, cacc];
            etimes{di, ci} = [etimes{di, ci}, etime];
            fvals{di, ci} = [fvals{di, ci}, fval];
            postdim{di, ci} = [postdim{di, ci}, K];
            macc{di, ci} = [macc{di, ci}, cmacc];
            mfval{di, ci} = [mfval{di, ci}, cmfval];
            aucs{di, ci} = [aucs{di, ci}, cauc];
            mauc{di, ci} = [mauc{di, ci}, cmauc];
            
            fprintf('(%d,%d,%d): K = %d, accu = %.4f, auc = %.4f, fval = %.4f, macc = (%.4f,%.4f), mauc = (%.4f,%.4f), mfval = %.4f\n', ...
                di, ci, loopi, K, cacc, cauc, fval, cmacc(1), cmacc(2), cmauc(1), cmauc(2), cmfval);
            
            if cachiter
                parsave([savedir filesep sprintf('opts_%d_%d.mat', di, ci)], ...
                    mZ, mW, meta, minvlambda, ...
                    accu{di, ci}, etimes{di, ci}, fvals{di, ci}, postdim{di, ci}, macc{di, ci}, mfval{di, ci}, aucs{di, ci}, mauc{di, ci}, ...
                    regvals, ell, alphav, invsigmaetasqr, invsigmawsqr, invsigmaxsqr, ...
                    poisstrunc, algtype, maxiter, burnin, max(cmacc) > maxacc);
            end
            maxacc = max(max(cmacc), maxacc);
            
            if algtype == 2
                if fval > minfval
                    cnvg = cnvg + 2;
                else
                    if all(deviation({minfval}, {fval}) < 1e-5)
                        cnvg = cnvg + 1;
                    else
                        cnvg = 0;
                    end
                    minfval = fval;
                end
            else
                if all(deviation({minfval}, {cmfval}) < 1e-5)
                    cnvg = cnvg + 1;
                else
                    cnvg = 0;
                end
                minfval = min(minfval, cmfval);
            end
            if loopi > 4 && max(cmacc) < accu{di, ci}(1,end-4)
                earlystop = earlystop + 2;
            elseif max(cmacc) < maxacc
                earlystop = earlystop + 1;
            else
                earlystop = 0;
            end
            if cnvg == 10 || loopi == maxiter || earlystop >= 10 || K > 2000
                break;
            end
            loopi = loopi + 1;
        end
        if ~cachiter
            parsave([savedir filesep sprintf('opts_%d_%d.mat', di, ci)], ...
                    mZ, mW, meta, minvlambda, ...
                    accu{di, ci}, etimes{di, ci}, fvals{di, ci}, postdim{di, ci}, macc{di, ci}, mfval{di, ci}, aucs{di, ci}, mauc{di, ci}, ...
                    regvals, ell, alphav, invsigmaetasqr, invsigmawsqr, invsigmaxsqr, ...
                    poisstrunc, algtype, maxiter, burnin, max(cmacc) > maxacc);
        end
    end
    save([savedir filesep 'opts.mat'], 'accu', 'etimes', 'fvals', 'postdim', 'macc', 'mfval', 'aucs', 'mauc', ...
        'regvals', 'ell', 'alphav', 'invsigmaetasqr', 'invsigmawsqr', 'invsigmaxsqr', ...
        'poisstrunc', 'algtype', 'maxiter', 'burnin');
end
